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Abstract. We have developed a three-dimensional radiative transfer method designed specifically for use with parallel adaptive 
mesh refinement hydrodynamics codes. This new algorithm, which we call hybrid characteristics, introduces a novel form of 
ray tracing that can neither be classified as long, nor as short characteristics, but which applies the underlying principles, i.e. 
efficient execution through interpolation and parallelizability, of both. 

Primary applications of the hybrid characteristics method are radiation hydrodynamics problems that take into account the 
effects of photoionization and heating due to point sources of radiation. The method is implemented in the hydrodynamics 
package FLASH. The ionization, heating, and cooling processes are modelled using the DORIC ionization package. Upon 
comparison with the long characteristics method, we find that our method calculates the column density with a similarly high 
accuracy and produces sharp and well defined shadows. We show the quality of the new algorithm in an application to the 
photoevaporation of multiple over-dense clumps. 

We present several test problems demonstrating the feasibility of our method for performing high resolution three-dimensional 
radiation hydrodynamics calculations that span a large range of scales. Initial performance tests show that the ray tracing part 
of our method takes less time to execute than other parts of the calculation (e.g. hydrodynamics and adaptive mesh refinement), 
and that a high degree of efficiency is obtained in parallel execution. Although the hybrid characteristics method is developed 
for problems involving photoionization due to point sources, the algorithm can be easily adapted to the case of more general 
radiation fields. 
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of our application of ionization calculations, the optical depth is 
used to determine the photoionization and heating rates. When 



Since astrophysical applications are many times dominated 
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fled s cenarios of explosions of massive stars ( Janka & Mueller 
1996), to name just a few. 

In creating a method that combines radiative transfer 
and hydrodynamics, one in general starts with an exist- 
ing hydrodynamics c ode and adds the ne c essary radiation 
proce ss es to it (e.g. [M ellema et al l Il998t I Turner & Stone 
| 200ll IWhitehouse & Bate! 120041 iHeinem ann et alJ 1 20051: 
iLiebendorfer et al .1120051) . In this paper we describe the addi- 
tion of a new radiative transfer algorithm, which we call hybrid 
characteristi cs, to the parallel 3 D AMR hydrodynamics pack- 
age FLASH (iFrvxell et alfeOOOh . 

Most of the radiative transfer methods that were success- 
fully combined with extant hydrodynamics codes apply some 
form of ray tracing to find the optical depth at each location 
in the computational domain. Apart from ray tracing one could 
also use statistical methods to find the solution to the radia- 
tive transfer equations (e.g. Ma selli etal1l2003h . Yet anothe r 
approach could be th e use of Fourier transf orms (lCenl l2002). 
or unstructured grids feitzerveld et al.lEo04l) . to determine the 
radiation field. 

For a one-dimensional, non-AMR, serial code, ray tracing 
becomes a rather straightforward procedure which requires lit- 
tle second thought. Equivalently, the case of a plane parallel 
radiation field on a Cartesian grid, or a single point- source at 
the centre of a spherically symmetric grid, for which all rays 
run parallel to a coordinate axis, can be handled quite easily. 
Although this type of implementation can readily be used to 
study a number of interesting astrophysical phenomena, it is 
still highly desirable to have a code that can treat the more 
general case of a point source of ionizing radiation on a 3D 
Cartesian domain. Such more general me t hods were for ex- 
ample implemented by Raeaet 

aD fcQQCt : iRichling & Y orke 
(2000); iLim & Mellemal (120031) . but none of these methods 
was explicitly parallelized for distributed memory machines 
though. 

The aim of this work is to create a characteristics-based ra- 
diative transfer method that can handle multiple sources of ion- 
izing radiation in AMR enabled simulations to be run on dis- 
tributed memory parallel machines. For this, a radical rethink 
of the concept of ray tracing is necessary, since, for this type 
of parallel AMR codes, the computational domain is not only 
sub-divided into a hierarchy of patches, but is also distributed 
over a number of processors. The first choice one therefore has 
to make is which flavour of ray tracing one wants to apply: ei- 
ther long or short characteristics. Since these two methods have 
rather different properties when it comes to efficiency and par- 
allelizability, this choice will determine the success of the final 
algorithm. 

We are aware of a number of other methods that use some 
form of adaptiv i ty to solve the radiative transfer equations: 
lAbel & Wandelt (2002) designed a method where the ray it- 
self is adaptive ly split into sub-rays, b ut the underlying grid 
is still regular. Steinacker et al. (2002) employed second or- 
der finite differencing of the full radiative t ransfer equations 
on an oct-tree AMR grid, and, more recently. I Juvela & Padoanl 
(120051) implemente d a ray tracing method for cell-based AMR. 
Ijessee et all Jl998h presented a radiative transfer method for 
patch-based AMR that uses the discrete ordinates approach. 




a b 



Fig. 1. Comparing the long (a) and short (b) characteristics 
method. For the long characteristics method, the closer one gets 
to the source, the more rays pass through (approximately) the 
same part of a cell, resulting in a large number of redundant 
calculations. The short characteristics method does not suffer 
from this, since here column densities are interpolated from 
cells that have been dealt with previously, so only the contribu- 
tions to the column density of the short ray sections that pass 
from cell to cell need to be computed. 

However, none of these methods resulted in parallel algorithms 
used in applications in which radiation is coupled to hydrody- 
namics. 

Efforts to create a parallel radiation hydr o dynam - 
ics code were p r esented by e.g. |Nakamoto et al. (2001); 
lHaves & Normanl (120031). and, more r e cently, a three- 
dimensional method by Heinemann et al. (2005), who de- 
veloped a ray tracing algorithm for decomposed domains. 
However, none of these two methods uses AMR. 

Our presentation begins with Sect. |2] in which we describe 
our new method. This method can not be classified as either 
short or long characteristics, but does have the desired prop- 
erties, namely high parallel and computational efficiency, of 
its two predecessors. We also compare our method to two re- 
cent ones which share similar features with ours. Supplemental 
physics components required by our primary target application 
(gas ionization, heating, and cooling) are presented in Sect. 
El where we give a brief de scription of the DORIC routines 
( Mell ema & Lundq vistll2002t and references therein). In Sect. 
|4| we first compare the accuracy with which our method cal- 
culates column densities to results obtained with a standard 
long characteristics approach. Then we present a pure radiation 
transport problem aimed at testing the accuracy of the ioniza- 
tion state calculations and shadow casting. This is followed by 
a coupled radiation hydrodynamics calculation of a photoevap- 
oration flow. Sect.|5]presents some initial performance results. 
Discussion of possible extensions and future applications for 
our method together with the conclusions are given in Sect.|6l 

2. Characteristics based radiative transfer 

When calculating the effects of ionizing radiation due to a point 
source, the radiation field is often dominated by this source, and 
one can safely ignore contributions to the radiation field due to 
the ambient gas. This means that the radiative transfer equa- 
tions assume a particularly simple form, since we can take the 
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total emission coefficient (and thereby the source function) to 
be equal to zero. Furthermore, when we also ignore the effects 
of scattering, the solution to the radiative transfer equations for 
the specific intensity I at location r is given by 

/(r) = /(0)exp(-r(r)), (1) 

and only depends on the optical depth r, which is defined by 

r(r) = a N(r) , (2) 

with ao the absorption cross section, and N the column density 
at r. 

Once the optical depth is known at every location in the 
computational domain, one can use it to find the ionization, 
heating, and cooling rates, and calculate the ionization state and 
temperature of the gas. Since, for finite- volume hydrodynam- 
ics codes, the computational domain is discretized into cells, 
the optical depth, or, equivalently, the column density for a cer- 
tain cell, is found by adding the contributions from all cells that 
lie between the source and the destination cell under consider- 
ation. This can be achieved by casting a ray, or long charac- 
teristic, from the source to the cell, accumulating contributions 
to the total column density along the way. In case of an AMR 
hierarchy, the algorithm first needs to identify the patches and 
cells contained within the patches that are traversed by the ray, 
and then calculate their local contributions to the total column 
density. 

Although the method of long characteristics is very ac- 
curate, it is also rather inefficient, since, the closer a cell is 
to the source, the more rays cut through (approximately) the 
same part of the cell, introducing a lot of redundant calcula- 
tions (see Fig. [l^). A way to eliminate this redundancy is to 
use the method of short characteristics. Here, the total column 
density for a certain cell is calculated by interpolating upwind 
values of column density calculated in a previous step, thereby 
creating some diffusion, but removing the redundant calcula- 
tions inherent in the long characteristics method (Fig.[H)). For 
this to work, the appropriate information from upwind cells 
needs to be available at all times, which means one needs to 
sweep the numerical grid outwards from the source. This ne- 
cessity of having to traverse the grid in a certain order makes 
this method intrinsically serial, since values of column density 
in cells now depend on one another. The long characteristics 
method does not suffer from this restriction, because here con- 
tributions to the total column density from cells cut by a ray do 
not depend on column densities in other cells. Therefore, the 
long characteristics method is fully parallelizable, since calcu- 
lations of contributions to the column density along each ray 
can be performed independently. For our method we combine 
the desirable qualities of both these approaches; the idea of in- 
terpolation is adopted from the short characteristics method, 
while parallelism is obtained following principles of the long 
characteristics method. 

In what follows, we start with a general description of the 
algorithm used to trace rays on AMR hierarchies. We explain 
how the long characteristics method is exploited to make this 
a parallel algorithm, and where the interpolation comes in to 
increase the efficiency of the calculation. 



processor processor 1 

Fig. 2. Two-dimensional example of an AMR hierarchy dis- 
tributed over two different processors. Here, each patch con- 
tains 4x4 cells. 

Although our algorithm is designed for three dimensions, 
many features of its implementation can be explained using 
two-dimensional analogues. Whenever the generalization from 
two to three dimension is non-trivial, we will supply the full, 
three-dimensional, description. Since the algorithm is naturally 
subdivided into a number of steps, we will expand on these sep- 
arately. 

2. 1. The distributed computational domain 

Consider a computational domain that is distributed over N p 
processors (for a two-dimensional example, see Fig.0. Rays 
are traced over these different sub-domains and must therefore 
be split up into independent ray sections. Naturally, these sec- 
tions are in the first place defined by the boundaries of each 
processor's sub-domain, and in the second place by the bound- 
aries of the patches contained within that sub-domain. 

So first each processor calculates for all the patches it owns 
the local column densities AN. These local contributions are 
found by tracing ray sections that originate at the patch faces 
that are located closest to the source, and that terminate at the 
centres of the cells (see Fig. 0. Since finding these contribu- 
tions is a local process, this part of the algorithm is fully par- 
allel, and can be implemented using either the short or long 
characteristics method. Details on how the ray tracing for in- 
dividual patches is implemented are given in Sect. I2.2I Note 
however that, before each processor can calculate its AN, it 
needs to know the physical location of the source, so this infor- 
mation is made available first. 

Since in general rays traverse more than one processor do- 
main, exchange of information has to take place at some point 
in the algorithm. After this communication step has finished, 
each processor should have available all contributions of col- 
umn density to the rays that terminate in its domain. By inter- 
polating and accumulating all these contributions for all rays, 
one ultimately obtains the total column density for each cell 
(see Sect. l2.3.Tl for a more elaborate description of the commu- 
nication patterns involved). Details on the procedure applied to 
find the patches cut by a ray, and the way in which their contri- 
butions to the total column density are subsequently calculated, 
are given in Sect. l2.3.2l and Sect. l2.3.3l respectively. 
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Fig. 3. Two-dimensional example of ray sections for a single 
patch. Local contributions to the column density are indicated 
by ray sections that terminate at cell centres (a), whereas con- 
tributions that are to be communicated between processors, and 
are subsequently used in an interpolation step, terminate at cell 
corners (b). The source lies outside of the patch in the direction 
of the lower left corner. 

2.2. Ray tracing a single patch 

In this section we will explain how the contribution to the total 
column density along a local ray section in a single patch can 
be calculated (see Fig. [3] and Fig.0). As explained above, each 
patch can be dealt with independently, which makes this part 
of the calculation fully parallel. 

The local column density contributions are calculated from 

A7V = ^x(HI)n(H)As, (3) 

cells 

with x(Bl) the ionization fraction of neutral hydrogen, n(H) 
the hydrogen number density, and As the physical path length 
through the cell. 

These contributions are found by casting a ray section from 
the faces of the patch that are located closest to the source to- 
wards each cell centre (see Fig. [3]). Column density contribu- 
tions by the cells that lie inside the patch along each section 
are calculated using the 'f ast voxel traversal algorithm' from 
lAmanatides & Wool Jl9871) (for more details on this traversal 
method, see App.lXl). Besides ray sections that terminate at cell 
centres, we also need to calculate the column density contribu- 
tion for ray sections that lead to cell corners located at those 
patch faces that are farthest away from the source (see Fig.[3j). 
These are the contributions to the column density that need to 
be communicated (Sect. 12. 3. Tl . and interpolated (Sect. 12. 3. "3b in 
subsequent steps of the algorithm. 

Calculating these ray sections is similar to the method of 
long characteristics, but since the number of cells per patch is 
low relative to the effective resolution of the full computational 
domain, this actually does not impair the performance of the 
method too much (see Sect. |3 for an analytical comparison of 
our method with the short and long characteristics one for a 
regular grid). 

We also considered using short instead of long characteris- 
tics to ray trace a single patch (see App. [0 for a description 
of a possible implementation). However, although the short 
characteristics method executes presumably more efficiently 



than the long characteristics one, the first requires interpola- 
tion, whereas the latter simply adds up column density con- 
tributions by individual cells. When the number of cells that 
need to be traversed is relatively small, as is the case when 
ray tracing the single patches, these extra calculations may ren- 
der the short characteristics method even less efficient than the 
long characteristics one. Furthermore, the interpolation intro- 
duces undesirable diffusion. We therefore decided to imple- 
ment the more accurate and straightforward ray tracing ap- 
proach of Amanatid es & Wool (Il987l) . 

2.3. Hybrid Characteristics 

As was mentioned above, in AMR hydrodynamics codes, each 
processor owns a sub-domain of the computational volume 
which is covered by a collection of patches. In order to obtain 
the total column density for a certain ray that traverses these 
sub-domains, individual local contributions by the patches need 
to be accumulated. This can be interpreted as applying the 
method of long characteristics, in this case not to add up con- 
tributions from individual cells, but instead to add up contribu- 
tions from individual patches. So here our algorithm does again 
make use of long characteristics but now at the level of patches. 

Since each processor knows the direction of its rays and the 
co-ordinates where they terminate, it can find the patches cut 
by these rays and perform the required calculations. For cer- 
tain flavours of AMR, patches from different refinement levels 
may partially overlap. In such cases, one would have to make 
sure that only parts of the patches that contain valid data (i.e., 
the data from regions resolved to the highest resolution) are 
considered in the calculation of the column density. One way 
to eliminate the overlap is to appl y a procedure calle d 'grid ho- 
mogenization', as described by Kre vlos et alJ J2002I) . 

For the oct-tree type of AMR implemented in FLASH, 
patches from different refinement levels do not overlap. Patches 
are either fully covered by still more refined patches or oth- 
erwise contain valid data (the latter are the so-called 'leaf 
patches' in terminology of FLASH). Therefore, a simple check 
to see if a patch is a 'leaf patch' is sufficient to determine 
whether or not it should contribute to the total column density 
along the ray. 

Once the list of patches traversed by a ray is known, we 
loop through it, and determine the local column density con- 
tributed by each patch to the total column density for the ray. 
Unless the ray terminates in the patch under consideration, it 
will in general not exit a patch exactly at a cell corner. This 
means that we need to interpolate the values of column den- 
sity contribution A7V, obtained earlier (using either the short 
or long characteristics method as described in Sect. 12.21) at that 
face of the patch where the ray leaves it. 

We would like to emphasize that, although our method 
makes use of some form of long characteristics, nowhere in 
the algorithm is a ray traced on a cell-by-cell basis over the full 
computational domain. To the contrary, ray sections are traced 
through the cells of each patch and it is these local contribu- 
tions which are combined through interpolation by performing 
another ray trace, this time not over cells but over patches, as 
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described in Sect. 12.3.21 This is why we call our algorithm hy- 
brid characteristics. 

Below, we first explain how the local column density con- 
tributions AN, obtained with one of the methods from Sect. 
12.21 are communicated between processors. Then we describe 
how the list of patches traversed by a ray is constructed, after 
which we show the way in which this list is used to calculate 
the contributions to the total column density N. 

2.3.1. Communicating local column density 
contributions 

Since, for a parallel AMR hydrodynamics code, the patches are 
distributed over a number of processors, communication be- 
tween processors is inevitable at certain points in the algorithm. 
In particular, as soon as the local contributions to the column 
density have been calculated (Sect. I2.2L values of these AN 
located at patch faces that are farthest away from the source are 
communicated between processors. In this way, each processor 
has the information regarding the face values of local column 
density from all patches in existence (i.e. the so called 'gather' 
operation is used). Apart from these face values, all processors 
also need information about the location and size of each patch 
and its refinement level in order to determine if a particular 
ray cuts a patch. This information is communicated using the 
'gather' operation as well. 

The size of the messages to be communicated and the mem- 
ory needed for storage of this information is given by 



processor 



processor 1 



-ftot Pmax S , 



(4) 



where P tot is the total number of processors, p ma x is the max- 
imum number of patches in existence on any processor, and S 
is the required storage space per patch. In three dimensions, S 
should contain the values of AN from the three patch faces 
located farthest away from the source, as well as the location, 
size, and refinement level information of each patch. 

For an initial test of the performance of the algorithm as 
a whole, and of its communication patterns in particular, see 
Sect.|3 

2.3.2. Constructing the list of patches cut by a ray 

A straightforward approach to constructing the list of patches 
traversed by a ray would be to simply check for all patches 
whether or not they are cut by the ray under consideration. 
Since this would have to be done for all rays, and since there 
are as many rays as there are cells, this approach quickly be- 
comes prohibitively slow. We therefore developed a new, more 
elaborate, but much faster method to find the list of patches cut 
by a ray. 

First, each processor creates a so called 'patch-mapping' 
which consists of an integer array representing the full compu- 
tational domain that stores the id (i.e. a unique integer identi- 
fier) of all patches containing valid data. In Fig.|4|we show an 
example of such a mapping. These local patch-mapping arrays 
then need to be communicated and merged (using a so called 
'reduce' communication operation) after which each processor 
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Fig. 4. Two-dimensional example of a 'patch-mapping' for a 
computational domain that is split over two different proces- 
sors. In the top row the local ids of the patches on the differ- 
ent processors are shown. The mapping of these patch ids onto 
the patch-mapping array is shown in the middle row. The bot- 
tom row shows the global patch-mapping after the local patch- 
mappings have been communicated. Tracing the depicted ray 
results in the patch list {1,4, 10, 13, 16}. The patch-mapping 
entries visited during the ray tracing are shown in grey. 

has the same global patch-mapping corresponding to the full 
computational domain. 

In order to discern patches that are on different processors 
we use the following coding for the global patch id: 



PG = PL +Ppn 



(5) 



with pc the global patch id, pl the local patch id, and P the 
processor id. 

We the n trace the ray, again using the 'fast voxel traversal 
algorithm' (Amanatides& see App.lAl. but now 

to trace through the global patch-mapping array. This results 
in the list of patches cut by the ray, which is used to accu- 
mulate their local contributions, which were already communi- 
cated earlier, to arrive at the total column density (as described 
in Sect. ESS . 

Although this approach to ray tracing can be a potential 
bottle-neck in the algorithm, one needs to keep in mind that 
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Fig. 5. Two-dimensional illustration of the linear interpolation 
scheme used to accumulate local column density contributions. 
Shown are the ray sections used in the interpolation (see text for 
further details). 



the maximum number of patch-mapping entries along a ray is 
given by V3C/c, with C 3 the total number of cells if the com- 
putational domain would be fully refined, and c 3 the number of 
cells per patch. 

For a typical three-dimensional oct-tree type AMR simula- 
tion with C = 512 and c = 16, we find a maximum amount of 
~ 55 patch-mapping entries that are cut by a ray. Note however 
that this is an upper limit. The number of entries is drastically 
smaller when the source and destination of the ray are not lo- 
cated at opposite sides of the domain (which will be the case for 
most rays). Note also that, although we have to trace through 
the patch-mapping entries, the actual number of patches that 
ends up in the list is strongly reduced due to the adaptive na- 
ture of the discretization. In the example given in Fig. |4] the 
number of patch-mapping entries visited by the ray is 13, but 
the number of patches that end up in the list is only 5. It is this 
latter number which determines how many interpolations are 
needed when accumulating the local column density contribu- 
tions. 

2.3.3. Accumulating local column density 
contributions 

Now that we have the list of patches traversed by a ray (Sect. 
12.3.21) and the values of local column density at the patch faces 
located farthest away from the source have been made available 
to all processors (Sect. 12. 3. Tl . we can proceed and calculate the 
local contributions to the total column density through interpo- 
lation. 

The calculations that need to be performed for a ray r 
traversing a patch p can be broken up into the following steps 
(two-dimensional case, see Fig.|3): 

1 Find the location e where r exits p. 

2 Use this to find the two cell corners c\ and c 2 that are clos- 
est to e and store their corresponding local column density 
contributions ATVi and A N 2 . 

3 Calculate the geometrical path lengths of the ray sections 
that terminate in ci, c 2 , and e, and denote these by li, l 2 , 
and l e , respectively. 



Fig. 6. Illustration of the interpolation scheme in three dimen- 
sions. For clarity we show outlines of cells on patch faces only. 
In the left image we show a ray r that exits the patch at location 
e through a cell face, together with the ray sections used in the 
interpolation that terminate at the corners of this cell face. The 
image on the right shows the cell face in more detail, where 
we indicated the cell corners by 1, 2, 3, and 4. In addition to 
these cell corners, ray sections used in the interpolation that 
terminate at 5, 6, 7, and 8 are also indicated (see text for further 
details). 

4 Use these path lengths to calculate the following normal- 
ized interpolation weights: 

Wt = \l 2 - l e \/(h + h), W 2 = \h - l e \/(h + l 2 ) . (6) 

5 Calculate the desired value of local column density at e 
through linear interpolation: 

AN e = w! ATVi + w 2 AN 2 . (7) 

After all A7V e for each patch in the list of patches cut by r are 
calculated, we simply need to sum them to arrive at the total 
column density for r: 

N(r) =Y, AN e(p) fe>elist(r)]. (8) 

v 

The interpolation weights given above were constructed us- 
ing the conditions 

wi li + w 2 l 2 = l e , and w\ + w 2 = 1 , (9) 

which, for the case of a homogeneous density distribution, re- 
sults in the exact solution for the column density (i.e., apart 
from a constant factor, the path length itself). Other weights, 
like ones derived from the distances between the exit locations 
e, ci, and c 2 , can also be used, but this leads to ~ 10% errors 
for rays that enter a patch close to a patch corner (as is depicted 
by the example ray section of Fig.|3). 

In three dimensions (see Fig.|6|) it is not straightforward to 
derive weights that are a generalization of the two-dimensional 
ones described above. We therefore give a more intuitive 
derivation of these weights, using a procedure where we apply 
the weights for the two-dimensional case twice in succession: 

1 Find the location e where r exits p. 
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Fig. 7. Summary of steps taken in the hybrid characteristics method. On the left we show ray sections that represent local 
contributions to the column density (summary step 3), whereas ray sections that represent values of column density that need to 
be communicated are shown at the centre image (summary step 4). Note that only those values on patch faces located farthest 
away from the source need to be communicated. On the right we show an example of the interpolation of these local values for 
a particular destination cell (summary step 6). Note that there is no need to interpolate the value for the final ray section in the 
destination patch since its value was already calculated previously (summary step 3). 



2 Use this to find the four cell corners ci, c 2 , cs, and C4 that 
are closest to e and store their corresponding local column 
density contributions AJVi, A7V 2 , A7V 3 , and A7V 4 . 

3 Calculate the geometrical path lengths of the ray sections 
that terminate in ci, c 2 , C3, C4, and e and denote these by 
h> h, h, h, and l e , respectively. Also calculate the path 
lengths h and Iq of the ray sections that terminate in C5 and 
cq respectively (see Fig.|6l>. 

4 Use these path lengths to calculate the following normal- 
ized interpolation weights: 

W! = \h ~ h\/(h + Z 2 ), w 2 = \h - h\/(h + l 2 ) , 
w 3 = |/4 - h\/(h + U), = \h - h\/(h + h) , 
w$ = \k - h\/(h + k), ^6 = I/5 - h\/(h + h) • 

(10) 

5 Calculate the values of local column density A A/5 and A A^ 
at C5 and cq through linear interpolation: 

A A/5 = Wl ANi + w 2 AN 2 , 

ANq = w 3 AN 3 + w 4 AA^ 4 . U } 

6 Calculate the desired value of local column density at e 
through linear interpolation of A A^5 and A A^ : 

AN e = w b AN 5 + w 6 AA^ 6 . (12) 

Our choice of using the values of local column density at 
C5 and cq to arrive at AA^ is arbitrary. Instead, one may also 



use the ones from cj and eg (cf. Fig.|6|) in the steps described 
above. 

The main difficulty in finding an interpolation scheme for 
the three-dimensional case lies in the fact that we need to weigh 
with the lengths of the ray sections to avoid the errors which 
will otherwise occur when the ray under consideration enters 
the patch close to a patch corner. Since in general all these path 
lengths are different from one another, this introduces quite 
a number of independent variables into the equations. So, al- 
though the two-step procedure just described is not unique, it 
is simple and fast, and it gives good results in practice. 

2.4. Summary of the algorithm 

The steps taken in the algorithm can be summarized as follows 
(see Fig.Q: 

1 Each processor checks if its sub-domain contains the 
source. The processor that owns the source stores its patch 
and processor id and makes it available to all other pro- 
cessors (broadcast). Note that this id may change during a 
simulation due to changes in refinement and the consequent 
redistribition of patches among processors. 

2 On each processor, create the local patch-mapping and 
communicate (reduce) it so that each processor ends up 
with the global patch-mapping (Sect. 12. 3. "2t . 
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3 On each processor, calculate local column density contri- 
butions AA^ for each patch using the 'fast voxel traversal 
algorithm' (Sect.Eland Fig.Qleft). 

4 Communicate (gather) all AN values at patch faces lo- 
cated farthest away from the source (Fig. centre). Also, 
the coordinates and refinement levels of all patches need to 
be gathered. This communication is done most efficiently 
when this information is combined in a single data type of 
size S (Sect. ETA . 

5 On each processor, construct for each ray the list of patches 
that are traversed by that ray f Sect. 12.3. "51 . 

6 On each processor, interpolate and accumulate the local 
contributions AN from the patches that are in the list to ar- 
rive at the total column density N (cf. Sect. 12.3.31 and Fig. 
□right) . 

2.5. Comparison to other methods 

To conclude this section, we compare our method to two more 
recent ones that either u se some form of adaptivity to trace rays 
(ljuvela & Padoanl2QQ5l), or that are parallelized for distributed 
memory architectures ( Heinema nn et alJ [20051) . Unlike ours, 
these methods are intended to solve for the full radiation field, 
and therefore need to employ multiple sets of rays to sample 
the angular parameter space. Depending on the adopted form 
of the source function, (lambda- iteration is to be performed as 
well in order to obtain a converging solution. 

ljuvela & Padoanl (120051) proposed a ray tracing method for 
cell based AMR intended to be used in calculations of line 
emission. Their method uses sets of parallel (in the geometri- 
cal sense) long characteristics to find the intensity at cell faces, 
which are then interpolated to get the intensity at the cell cen- 
tre using a short characteristic. This is repeated for a number of 
directions after which angle averaged quantities are obtained. 
This process is then lambda-iterated to get converging line in- 
tensities. 

Since their method refines on a cell-by-cell basis, and ours 
employs patches structured in an oct-tree hierarchy, there is a 
one-to-one correspondence between the procedures of ray trac- 
ing used in the two methods: their long characteristics corre- 
spond to our ray tracing of the patch-mapping, whereas their 
short characteristics correspond to our ray tracing of a single 
patch. 

More recently, Heineman n et al ] §005) developed a 
method for tracing rays through a decomposed computational 
domain (i.e. sub-domains that are distributed over a number of 
processors). To sample the radiation field, rays are traced that 
are either parallel or diagonal to a regular patch. As they men- 
tion, this means that there is no need for them to interpolate 
local values. Furthermore, since their source function acts only 
locally, their is no need to iterate the solution. 

As in our approach, Hei nemann et alJ (120051) first obtain all 
local contributions (which they call 'intrinsic') and add these 
up to arrive at the total solution. Ho wever, in contrast to our 
method, the communication pattern of Heinem ann et alJ (120051) 
is intrinsically serial (i.e. processors have to wait for one an- 
other, see their Fig. 1). In their case of a decomposed regu- 



lar domain, the performance penalty due to the serial nature of 
their algorithm is small, but in case of an A MR type of grid, 
the per formance would be severely degraded. Heineman n et al I 
( 2005) also consider the special case of periodic boundary con- 
ditions with rays running only parallel along a coordinate axis. 
In such a situation, the boundary values are broadcasted and 
the inter-processor communication is more efficient than in the 
serial case. 

Although our method is designed to study the effects of 
ionization due to point sources of radiation, it can be eas- 
ily adapted to trace sets of parallel rays instead. Depending 
on the application, a prescription for the source function and 
(lambda-)iteration would need to be implemented. This would 
make our method suitable for solving the radiative transfer 
equation in a more general way, similar to the methods just 
discussed. The added advantage of such an approach is that our 
method is highly parallel and coupled to an AMR hydrodynam- 
ics code. 



3. Ionization, heating, cooling 

When the column density from the source up to each cell face 
is known, the ionization fractions and temperature can be com- 
puted. For th is we use a simplified version of the D ORIC 
routines (see Mell ema & Lundqvistl l2002t iFrank & Mellema 
1994). In what follows, we summarize the way in which these 
routines calculate the ionization, heating, and cooling rat es 
(for more details, please refer to Frank & Mellema (1994)). 
Although the DORIC package is capable of handling a large 
number of species (H, He, C, N, O, and Ne), we use hydrogen 
as the only gas component in order to keep the complexity of 
our tests cases at a minimum, and we will therefore describe 
just this case. 

The ionization fractions of hydrogen are given by 



with 



*(hd - .(HID = 2m 

[ ' n(H) ' [ ' n(H) ' 



n(H) = n(HI) + n(HII) 



(13) 



(14) 



the total hydrogen number density. The electron number den- 
sity follows from 



n(HII) +n(C), 



(15) 



where the number density of carbon is included to prevent the 
possibility of n e = 0, by assuming that carbon is always at 
least singly ionized due to the interstellar UV field. 

For hyd rogen, the n umber of photoionizations per second 
is given by (Osterbrock 1989) 



ao dv, 



(16) 



with J v the local mean intensity of the radiation field, ao = 
6.3 x 10 -18 cm 2 the cross section (which we take to be fre- 
quency independent, or 'grey', for simplicity) and v§ the ion- 
ization threshold frequency. 
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The number of collisional ionizations per second is calcu- 
lated using 

A c = A c (m)n e VTexp(-I(m)/kT) (17) 

with A c (HI) = 5.84 x lO^cmfK- 1 / 2 , and /(HI) the hydro- 
gen ionization potential (Cox 1970). 

For the on-the-spot appr oximation, the radiative recombi- 
nation rate is given by (cf . lOsterbrockl 19891) 

-0.7 



a R = a R (10 4 K) 



io 4 7 



(18) 



13 cm 3 s 1 . The temperature is 



with a R (10 4 K) = 2.59 x 10 
determined from the pressure using 

p= (n(R)+n e )kT. (19) 

The rate equation for the hydrogen ionization fraction is 
given by 

dx(HII) 



dt 



c(m)A - x(HII)n e a R , 



with A = A p + A c the total number of photo- and collisional 
ionizations per second. 

When one assumes that the electron density is constant, an 
analytic solution for x(HII) can be found, and iterating for n P 
gives the time dependent solution ( Schmi dt- Voigt & Koeppen 
1987). Since A c and a R are both temperature dependent, the 
change in temperature due to heating and cooling needs to be 
recalculated for each iteration step as well. 

The photoionization heating rate is given by 



ra(HII) 



hv 



-a§h(y — z/q) dv, 



(21) 



and for the coolin g rate we use a coll i sional equilibrium 
cooling curve from Dalearno & McCravl (1 19721) (more gen- 
eral composition-dependent cooling is available in the DORIC 
package). 

For calculating the local mean intensity of the radiation 
field we use a blackbody spectrum, so we have 



47rJ„(r) 



hv 3 



■exp(-r(r)). (22) 



Here, Rs is the radius, and T$ is the effective temperature of 
the source. The optical depth r at position r of equation © is 
calculated using the total column density N(r) from equation 
(|8|). Since evaluating the integrals for the photoionization and 
heating rate [equations (fT6t and (l2Tl l is too time consuming to 
perform for every value of r, they are stored in look-up tables 
for a range of optical depths and interpolated when needed. 

The hydrodynamics and ionization calculations are cou- 
pled through operator splitting. To avoid having to take time 
steps that are the minimum of the hydrodynamics, ioniza- 
tion, and heating/cooling time scales, we use the fact that the 
equations for the ionization and heating/cooling can be iter- 
ated to c onvergence. Since these are so called 'stiff equa- 
tions (e.g. |Press_et alJI 19921) . we use a special iteration scheme 
( Frank & Mellema 1994|)7This means that the only restric- 
tion on the time step co mes from the hydrodynam ics (i.e. the 
Courant condition). See Frank & MellemjJ i 19941) for an as- 
sessment of the validity of this approach. 



4. Tests 

In this section we present a number of tests for our new al- 
gorithm. First, we discuss the accuracy with which column 
densities and ionization fractions are calculated. We compare 
the results obtained with the hybrid characteristics method to 
those calculated using long characteristics. Since the interpola- 
tion scheme of our method is designed to give the exact result 
for the column density in case of a homogeneous density dis- 
, we also consider its performance in case 



tribution fSect. 12.3.31) 

of a more general density field. 

We conclude this section by testing the shadow casting ca- 
pabilities of our method and apply it to a 'real-world' appli- 
cation of photoevaporating flows. This last test demonstrates 
the performance of out method when used in combination with 
hydrodynamics. 



(20) 4. 1. Column density 



We performed two-dimensional calculations where we placed a 
single point source at the centre of a 1 /r 2 density distribution, 
the result of which is shown in Fig. [8] In order to prevent an 
under-resolved singularity at the location of the source, we used 
a constant density sphere with a radius of 5 x 10 14 at the source 
location. 

The left panel of Fig. [8] shows the column density distri- 
bution along a line y = const cutting through the exact loca- 
tion of the source. Since for this special case no interpolation is 
necessary, only very small differences between the two meth- 
ods are found. These differences are due to uneven sampling of 
the 1/r 2 density distribution on the adaptive mesh. The errors 
increase only slightly (< 0.5%) when interpolation is used, as 
indicated by the results obtained along the y = const line lo- 
cated at 3/4 of the horizontal extent of the domain (the right 
panel in Fig. [8]). 

4.2. Shadow casting 

To test the shadow casting capabilities of our algorithm, we cal- 
culate the column density and ionization fractions for a homo- 
geneous environment containing higher density clumps, which 
are taken to be spherical and neutral. The ionization state is 
found by iterating the ionization fractions over a period equal 
to a few recombination time scales, while keeping the temper- 
ature fixed. 

The computational domain spans the region 
(2.0,1.0,1.0) x 10 18 cm. The ambient medium has a 
number density n env = 10 2 cm -3 and a temperature 
T env = 5000 K. The source of ionizing radiation is located 
at (x,y,z) = (0.0,0.5,0.5) x 10 18 cm. It has a luminosity 
L s = 7000 L© and an effective temperature T eff = 50000 K. 
The resulting Stromgren sphere has a radius that is larger 
(~ 3 x 10 18 cm) than the physical size of the computational 
domain. Two identical clumps are placed at a distance of 
~ 10 18 cm from the source. Each clump has a density 
^ciump = 10 4 cm -3 , a temperature T c i ump = 100 K, and a 
radius r c i ump = 4 x 10 16 cm. We used 6 levels of refinement 
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Fig. 8. Values of column density for the case of a single point source in a two-dimensional domain with a 1/r 2 density distribu- 
tion. Shown are one-dimensional cuts along the y-direction through the source located at the centre of the domain (left panels) and 
at 3/4 of the domain (right panels). In the top panels, the solid line indicates the result for the long, whereas the crosses indicate 
the result for the hybrid characteristics method. The bottom panels show the ratio (hybrid/long) of column density values. 



with patches of 16 3 cells. The effective resolution in this test 
was 1024 x 512 2 cells. 

The results of the shadow casting test are shown in Fig. |9] 
As one can see, our hybrid characteristics method is capable 
of casting shadows with very sharp boundaries, indicating a 
low numerical diffusivity of the scheme. We note that since the 
initial conditions do not contain any density gradient, column 
densities calculated in this test are identical to the ones one 
would obtain using a long characteristics method. 



4.3. Application: photoevaporating clumps 

To illustrate that our hybrid radiative transfer algorithm can 
be used efficiently in combination with hydrodynamics, we 
present a first 3D application of the evolution of over-dense 
clumps being photoevaporated. We use the parameters of the 
simulation setup described in Sect. l4.2l as initial conditions and 
follow the dynamical evolution for ~ 4000 yr. This simula - 
tion is similar to the ones presented by Lim &Mellema(2003), 
with this difference that in our simulation both the source and 
the clumps are inside the computational domain, and that our 
radiation field is not approximated by parallel rays. 

These computations are relevant to the shaping and evolu- 
tion of cometary knots which are observed in objects like the 
Helix (NGC 7293), Eskimo (NGC 2392), and Dumbbell (M27) 
nebula. Another application is the interaction zone that is ob- 
served between binary proplyds in HII regions like NGC 3603 



jBrandner et alJ Eoool) and the Orion Nebula ( Graha m et alJ 
l2002h . 

Fig -[El shows a sequence of snapshots 1 of the density and 
neutral hydrogen fraction at different times during the simu- 
lation. One sees that the interaction of the photoevaporation 
flows coming from the clumps results in a zone of higher 
density between the clu mps, which, as was already found by 
iLim & Mellema (2003), can explain the excess emission ob- 
served between some cometary knots and binary proplyds. This 
interaction zone recombines, becomes optically thick, and casts 
a shadow. It is interesting to see that the zone, and the shadow 
region behind it, persist even after the two clumps have been 
fully evaporated. This mechanism for creating extra shadows 
may influence the evolution and survival time of clumps that 
lie farther away from the star, an effect not taken into account 
in previous numerical studies. We are currently investigating 
further into this kind of flows and will present our findings in a 
forthcoming publication. 

5. Performance analysis 

We start this section by comparing our hybrid characteristics 
method to the more traditional long and short characteristics 
methods for regular grids. In order to do this, we distinguish 
between two types of computations: first we determine how 
many calculations are needed to arrive at the local contribution 



1 Movies of the simulations are available with the electronic version 
of this paper at http : / /www . edpsciences . org 
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Fig. 9. Values of log 10 of the column density (left) and log 10 of the HI ionization fraction (right) for the case of a single point 
source in an environment with a homogeneous density distribution containing neutral clumps with higher density. Shown are 
color coded plots of xy-cuts through the centre of the domain (top row) and xz-cuts through the centre of the bottom clump 
(bottom row). The bottom left image shows the AMR patch distribution superimposed, where each patch contains 16 3 cells. 



each cell makes to the column density along a ray, and second 
we look at the number of interpolations the different methods 
have to perform to compute the total column density up to each 
cell. 

Consider a computational domain with a resolution of C 3 
cells and a source located at one of the corners of the domain. 
For the case of a regular grid, the maximum number of cells 
a long characteristic would encounter is V3C, and, since we 
assume that a ray is cast to all cells, the number of calculations 
needed to provide the total column density is therefore < C 4 . 

For our hybrid characteristics method, which employs an 
oct-tree type of AMR grid, the maximum number of cells a lo- 
cal ray section encounters is y/3c, where c 3 is the number of 
cells in a single patch (cf. Sect. 12.21) . So in this case, for a fully 
refined grid, the total number of calculations would amount to 
< cC 3 . But, since in general the domain would be refined by 
a factor r, with < r < 1, this number reduces to < r cC 3 . 
This means that when rc~l, our method needs ~ C 3 calcu- 
lations to arrive at the local contributions to the column density 
for each cell, which is of the same order a short characteristics 
method would need on a regular grid. However, the local values 
of column density still need to be communicated and interpo- 
lated to arrive at the total column density for each cell. On the 
other hand, a short characteristics method also needs to inter- 
polate local values when it sweeps through the grid, whereas 
a long characteristics method, although it executes a factor C 
more calculations, does not need to perform any interpolations 
at all. 

The number of interpolations to be performed by our 
method is determined by the number of patches that are en- 
countered when ray tracing through the patch-mapping (cf. 
Sect. 12.3.21 and 12.33b . This number is at most V3C/c, since, 
for a fully refined grid, there are C/c patches along a coordi- 
nate axis. For a grid that is not fully refined this number is again 



reduced by a factor r. A ray trace through the patch-mapping 
is to be performed for every cell, which brings the total num- 
ber of interpolations to < r 2 (C / c)C 3 , where one factor of r 
comes from the number of patches cut by a ray, whereas the 
other factor comes from the total number of cells that exist in 
the computational domain. 

A short characteristics method needs to do an interpolation 
for every cell, so, for a regular grid, the total number of interpo- 
lations is C 3 . This implies that when r 2 C/c ~ 1, our method 
needs to compute a similar number of interpolations as a short 
characteristics one. Note that we assume that the calculations 
needed to do the interpolations are comparable in execution 
speed for the short and hybrid characteristics methods, which 
may actually not be the case. 

As an example, a typical AMR calculation has C = 512, 
c = 16 (i.e. 6 levels of refinement), and r = 0.25, which re- 
sults in rc — 4 and r 2 C/c — 2. This shows that, for a single 
processor calculation with a proper choice of the ratio C/c and 
a reasonable amount of refinement, our hybrid characteristics 
method is expected to perform equally well as a short char- 
acteristics method on a regular grid. It also means that, when 
our method is used in parallel, a better performance will be ob- 
tained when increasing the number of processors. 

To investigate this aspect in some more detail we have con- 
ducted a preliminary performance analysis using the photoe- 
vaporating clumps test case described in Sect. l4.3l as the under- 
lying physics problem. We used 5 levels of refinement irrespec- 
tively of the number of processor used in the test run (i.e., the 
problem had a fixed total work). Calculations have been ter- 
minated after reaching 10% of the nominal simulation time. 
Otherwise the simulation parameters were identical to those 
used in the calculations presented in Sect. 14.31 

The results of our performance study are shown in Fig.fTTl 
where we compare the overall performance of the hydrodynam- 
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Fig. 10. Snapshots of the evolution of the log 10 of the mass density (left) and the log 10 of the HI ionization fraction (right) for 
the case of a single point source in an environment with a homogeneous density distribution containing two neutral clumps with 
higher density. The source is located at (0.0, 0.5, 0.5) x 10 18 cm. Shown are color coded plots of xy-cuts through the centre of 
the domain at t = yr (first row), t = 792 yr (second row), t = 1584 yr (third row), t = 2377 yr (fourth row), t = 3169 yr (fifth 
row), and t = 3961 yr (sixth row). 
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Fig. 11. Performance of the main components of a radiation 
hydrodynamics calculation. 
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Fig. 12. Performance of the different steps in our hybrid char- 
acteristics method. For this specific test, the communication 
takes as much time as the rest of the calculation when using 
64 processors. As is explained in the text, for patches with a 
larger number of cells, this constraint may become less severe. 



ics, AMR, and radiation calculations. Detailed results for the 
radiation part are presented in Fig. [l2l where the timings for 
the individual components of our hybrid characteristics method 
(local ray trace, communication, accumulation/interpolation, 
and ionization) are shown. 

Performance data obtained for our realistic test problem in- 
dicates that the ray tracing part of the calculation takes less time 
than either the hydrodynamics or grid adaptation. Furthermore, 
it shows that most of the computational time during ray tracing 
is spent in interpolating and accumulating the local contribu- 
tions to the column density (i.e. step 6, Sect. 12.41) . Following the 
analytical assessment made above, we conclude that in order to 
reduce the number of interpolations required during calcula- 
tion one should try to minimize the value of r 2 C/c rather than 
rc when setting up a simulation. This suggest that one should 
use patches that contain a relatively large number of cells com- 
pared to the effective resolution of the computational domain, 



and, of course, keep the filling factor of the finest AMR level at 
a minimum. 

Fig-Elshows performance results for the radiation module 
including the ionization package for our fixed size problem. As 
one can see, the time required to calculate the column densities 
is about the same as the time needed to calculate the ionization 
state of the gas. Furthermore, both these calculations are local 
and therefore perform very well. On the other hand, we notice 
that the communication part of our algorithm does not perform 
perfectly. This is somewhat expected since, with the increasing 
number of processors, the efficiency of our algorithm becomes 
limited by the efficiency of the global gather operation (used to 
collect column densities from patch faces). The results for this 
specific test indicate that communication is likely to dominate 
the runtime when more than ~ 64 processors participate in the 
computations. We expect that this limitation becomes less se- 
vere when larger patches are used in the simulation. In this case 
the cost of communication may still be lower than, for exam- 
ple, the time needed to accumulate and interpolate the column 
densities. To determine whether this is indeed the case, more 
elaborate performance tests involving a larger number of pro- 
cessors (of the order of ~ 1000) are required, and such tests are 
currently underway. 

6. Conclusions 

We described a new radiative transfer algorithm for parallel 
AMR hydrodynamics codes, called hybrid characteristics. We 
presented details of several aspects of the algorithm: ray trac- 
ing, communication, and interpolation. 

The ray tracing is performed in two steps. First, local long 
characteristics are used to calculate column density contribu- 
tions for each patch. A second ray trace is then performed 
where a so called patch-mapping is used to find the patches cut 
by each ray. When the list of patches cut by a ray is known, in- 
terpolation of local column density values is required to find the 
total column density up to each cell. For this, one needs the val- 
ues of local column density contribution at patch faces, which 
are communicated to all processors. The coefficients used in 
the interpolation are chosen such that the exact solution for the 
column density is retrieved when there are no gradients in the 
density distribution. 

For the case where the distribution is not constant but has 
a 1/r 2 profile we find deviations of the order of ~ 0.5% when 
comparing our method with a long characteristics one. This 
high accuracy with which column densities values are calcu- 
lated results in well defined and sharp shadows. 

We showed that our method can be used efficiently for par- 
allel radiation hydrodynamics calculations in three dimensions 
on AMR grids. We presented preliminary results for our new 
method in application to the problem of the photoevaporation 
of two over-dense clumps due to the ionization by a single 
source of radiation. The results of this simulation offer a pos- 
sible explanation for the excess emission observed in between 
cometary knots seen in for example the Helix Nebula, and the 
interaction zone observed in binary proplyds found in HII re- 
gions like NGC 3603 and the Orion Nebula. These simulations 
also suggest a possible mechanism for the creation of extra 
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shadows by the high density interaction zone forming between 
the clumps. This additional shadowing may influence the evo- 
lution and survival time of clumps that lie farther away from 
the source. We are currently investigating further into this kind 
of interactions between photoevaporating flows, and their con- 
sequences for the dynamics, and will report our findings in a 
future publication. 

An initial performance test showed that our method works 
very well when used for calculations on a parallel machine. 
For this specific test, the communication part of our algorithm 
starts to dominate the calculation when more than ~ 64 proces- 
sors are used. However, we showed analytically that a careful 
choice of the ratio of the number of cells per patch to the to- 
tal number of cells in the computational domain controls the 
amount of communication used in the calculation. This anal- 
ysis can be used to optimize the design of our method. More 
in-depth performance and scaling studies are currently under- 
way, using large (~ 1000) number of processors, and these will 
also be used to further optimize the current implementation. 

Because of the modular nature of the FLASH code and the 
DORIC routines, additional elements like more sophisticated 
cooling or multiple species can easily be added. Also, multi- 
ple point sources can be handled by our method, and in prin- 
ciple moving sources could be implemented. Furthermore, it 
should be straightforward to extend the hybrid characteristics 
method so that it can be used to solve for a more general radia- 
tion field, with the added advantage that our method is already 
parallelized and coupled to an AMR hydrodynamics code. 

Another possible application for our method is the calcu- 
lation of the propagation of ionization fronts in a cosmolog- 
ical context. For these cal culations photon con servation is an 
important issue. Recently, iMellema et alJ (120051) developed a 
method for following R-type ionization fronts that may move 
more than one cell per time step, where a special formulation of 
the equations ensures photon conservation. Although the paral- 
lel nature of our algorithm may complicate the implementation 
of such an approach, we may s till benefit from the ideas pre- 
sented by Me llemaet alJ J2005h . 

We intend to make our method publicly available in a future 
FLASH release. 
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Appendix A: Fast voxel traversal 

Here we briefly discuss the 'fa st voxel traversal algorithm' 
fromlAmanatid es & Wool Jl987l> . We have used this algorithm 
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Fig. A.l. Explanation of different quantities used in the fast 
voxel traversal method. Shown is a single patch with a local 
long characteristic ray section. 

twice in our method, once to ray trace through the cells ('vox- 
els') of a single patch (local long characteristics, Sect. 12.21) . 
and once to ray trace the patch-mapping (hybrid characteris- 
tics, Sect. 12.31) . The idea behind the algorithm is to keep track 
of three different ray parameters t x , t y , and t Z9 one for each 
co-ordinate direction, and to use these to determine how to step 
from cell to cell through the patch, ensuring that all cells cut 
by the ray are visited (see Fig. IA.ll> . First, values for the incre- 
ments in ray parameter t needed to step from cell to cell in the 
x-, y-, and ^-direction, indicated by St x , St y , and St z , respec- 
tively, are determined: 

St x W/Ax ,St y t max /Ay ,St z t max /Az , (A.l) 

where t max = y Ax 2 + Ay 2 + Az 2 is the final ray parameter 
(i.e. the total path length of the ray). Next, the ray parameters 
t x ,t y , and t z are set to their respective initial values, indicated 
by t x ,ty, and t\, after which a loop is entered where the min- 
imum of these three values is determined. This gives the co- 
ordinate direction in which the cell lies that is to be visited next 
by the ray. For example, if min(^, t y ,t z ) = t x , the next cell 
the ray will enter lies in the x-direction, and we have to in- 
crement the ray parameter for the x-direction accordingly, i.e. 
t x = tx + St x . We loop as long as all ray parameters are smaller 
than the final ray parameter £ max . As a by-product, the algo- 
rithm produces the path length of the ray section for each cell 
that is crossed, which is obtained by subtracting the previous 
from the current ray parameter. 

Appendix B: 

Ray tracing a single patch: short characteristics 

As an alternative to the 'fast voxel traversal algorithm' for ray 
tracing a single patch as presented in Sect. 12.21 we here briefly 
describe the short characteristics method, which could be used 
for the same purpose. Since the method of short characteristics 
uses interpolation from neighbouring cells, upwind values need 
to be available at all times, so cells need to be swept in a certain 
order. This sweeping sequence is determined by the physical 
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Fig. B.l. Two examples of short characteristics sweeping se- 
quences for a single patch that may occur in practice. For the 
illustration on the left, the source is located inside the patch at 
the starting point of curve 1. In this case the four space filling 
curves shown should be swept in the order indicated. For the 
patch on the right, the source is external to the patch, and lies 
in the direction of the lower left corner, so there is only one 
curve that needs to be swept. 





Fig. B.2. Illustration of the interpolation scheme for a single 
cell used in the short characteristics method for the 2D (left) 
and 3D (right) case. 

location of the patch relative to the source position (see Fig. 
IB .1ft . Using the known physical location of the source, the ge- 
ometrical path length of the ray section that crosses a cell is 
calculated for every cell contained in the patch. The short char- 
acteristics method then sweeps the patch in a direction away 
from the source, interpolating upwind column density contri- 
butions for each cell along the way. 

For the two-dimensional case, Fig. IB. 21 illustrates which 
two cells, indicated by c\ and c 2 , are used in this interpolation. 
Simple linear interpolation weights 



w\ = 1 — d; w 2 = d 



(B.l) 



could be used to arrive at the column density contribution at 
cell c, using 



AN C 



HANi + Am . 



(B.2) 



with ANi the upwind values of column density that need to 
be interpolated, d the normalized distance from c\ to the loca- 
tion where the ray pierces the line connecting c\ and c 2 , At 
the physical path length of the short characteristic ray section, 
and n the number density inside the cell crossed by this short 
characteristic. 



For the three-dimensional case the ray pierces a cell face, 
so four instead of two quantities need to be interpolated. The 
normalized weights are chosen to correspond to the partial ar- 
eas of the cell face defined by the corners of this face and the 
location at which the ray leaves the cell (cf. Fig. IB. 2b : 

w x = (1 - di)(l - d 2 ); w 2 = di (1 - d 2 ); 
w 3 = (1 - di) d 2 ; = di d 2 . 



(B.3) 
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